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ALIAS-FREE  SMOOTHED  WIGNER  DISTRIBUTION 
FUNCTION  FOR  DISCRETE-TIME  SAMPLES 

INTRODUCTION 

The  utility  of  the  Wigner  distribution  function  ( WDF )  for 
detailed  time-frequency  analysis  of  waveforms  has  been  summarized 
very  well  in  a  recent  article  by  Cohen  [1];  this  material  will  be 
assumed  to  be  known  by  the  reader.  As  for  actual  numerical 
calculation,  the  problem  of  obtaining  an  alias-free  WDF  and 
complex  ambiguity  function  (CAF),  from  discrete-time  samples,  was 
solved  in  a  recent  report  by  Nuttall  [2].  Specifically,  an  upper 
bound  on  the  time  sampling  increment  and  a  lower  bound  on  the 
fast  Fourier  transform  ( FFT )  size  were  determined  that  allowed 
for  evaluation  of  the  original  continuous  WDF  and  CAF  at  a 
discrete  set  of  points  with  sufficient  detail  and  coverage  to 
avoid  any  significant  loss  of  information.  Furthermore,  a 
detailed  prescription  for  the  required  data  processing  of  the 
available  discrete-time  samples,  in  terms  of  FFTs ,  was  given. 

However,  the  presence  of  large  oscillating  interference 
terms,  which  are  inherent  to  the  WDF,  requires  that  some  smoothed 
version  of  the  WDF  be  made  available  from  discrete  data.  This 
problem  was  addressed  recently  by  Harms  [3],  and  a  procedure  was 
delineated  for  its  realization  in  terms  of  FFTs.  However,  the 
additional  data  processing  required  for  the  smoothed  WDF  cannot 
be  realized  without  some  extra  effort  or  penalty;  in  fact,  new 
more  stringent  bounds  on  the  sampling  increment  and  FFT  sizes 
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must  be  met  in  order  to  retain  the  alias-free  character  of  the 
resultant  smoothed  WDF.  These  bounds  were  derived  by  Nuttall  and 
furnished  to  Harms  who  listed  them  in  [3;  section  4  (see  his 
reference  11)]. 

In  this  current  report,  we  will  present  the  detailed 
derivations  that  lead  to  these  bounds.  In  the  process, 
interpretations  of  the  smoothed  temporal  correlation  function 
( TCF )  and  smoothed  spectral  correlation  function  (SCF)  are 
required  and  furnished.  Allowance  for  a  very  general  form  of 
ambiguity  weighting  (multiplication)  or  Wigner  smoothing 
(convolution),  including  tilts  in  the  appropriate  time-frequency 
planes,  is  made  and  accounted  for.  The  specific  data  processing 
and  FFT  operations  are  presented  in  complete  detail. 
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CONTINUOUS  TIME-FREQUENCY  REPRESENTATIONS 

In  this  section,  waveform  s(t)  is  considered  to  be  available 
for  continuous  time  t.  We  will  point  out  some  basic  properties 
of  the  various  time-frequency  representations  (TFRs)  of  the 
waveform,  which  will  be  required  later  when  we  address  the 
discrete-time  case;  some  of  this  material  was  given  in 
[2;  especially  appendix  A]. 


WAVEFORM  CHARACTERISTICS 

Complex  waveform  s(t)  has  voltage  density  spectrum 

S(f)  =  J  dt  exp(-i2nft)  s(t)  ,  (1) 

where  f  is  cyclic  frequency  and  integrals  without  limits  are 
conducted  over  the  range  of  nonzero  integrand.  It  will  be 
presumed  that  the  waveform  is  essentially  time  limited  and 
frequency  limited;  that  is, 

| s ( t ) |  *  0  for  | t |  >  T/2  (2) 

and 

| S ( f ) (  *  0  for  jf I  >  F/2  .  (3) 

Thus,  the  total  time  extent  of  s(t)  is  T  seconds  while  the  total 
frequency  extent  of  S(f)  is  F  Hertz.  The  effective  extent  of 
s(t),  say  where  |s(t)|  is  within  1/e  of  its  peak,  is  smaller  than 
T;  similarly,  the  effective  extent  of  S(f)  is  smaller  than  F. 
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This  distinction  between  the  essential  (total)  extent  and  the 
effective  extent  is  kept  below.  The  time-bandwidth  product  TF 
must  be  larger  than  1  and  can  be  much  larger  than  1  for  some 
waveforms  with  detailed  amplitude-  and/or  frequency-modulation. 

The  fact  that  s(t)  is  centered  at  t  =  0  results  in  no  loss  of 
generality  because  we  can  delay  or  advance  a  given  waveform  to 
this  location.  Similarly,  a  centered  spectrum  S(f)  is  easily 
achieved  by  frequency  shifting.  We  allow  for  complex  svc), 
thereby  accommodating  analytic  or  complex  envelope  waveforms. 


TIME-FREQUENCY  REPRESENTATIONS 

The  temporal  correlation  function  ( TCF )  of  s(t)  is  defined  as 

R  ( t ,  t  )  =  s  ( t+Jjr  )  s^t-^r)  .  (4) 

Reference  to  (2)  immediately  reveals  that  R(t,r)  is  essentially 
confined  to  |t|  <  T/2,  |t|  <  T.  The  quantity  r  is  the  time  delay 

or  separation  variable. 

The  spectral  correlation  function  (SCF)  is  the  double  Fourier 
transform  of  R(t,t)  and  is  given  by 


*( v, f  ) 


dt  dr  exp( - i2nvt-i2nfr  )  R(t,r)  = 
=  S(  f+*av )  S*(f-»sv)  . 


(5) 


Use  of  (3)  then  demonstrates  that  $(v,f)  is  essentially  limited 
to  |v|  <  F,  |f|  <  F/2.  The  quantity  v  is  the  frequency  shift  or 
separation  variable. 
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The  Wigner  distribution  function  (WDF)  is  then  given  by 
either  of  the  following  transforms 

W(t,f)  =  J  dt  exp(-i2nfx)  R(t,r)  =  (6a) 

=  J  dv  exp(+i2nvt)  #(v,f)  .  (6b) 

From  (6a),  we  can  conclude  that  W(t,f)  is  confined  to  |t|  <  T/2, 
while  from  (6b),  the  frequency  extent  is  essentially  |f|  <  F/2 . 

Finally,  the  complex  ambiguity  function  (CAF)  is  available 
from  either  of  the  following  transforms 

X(v,t)  =  J  dt  exp(-i2nvt)  R(t,r)  =  (7a) 

=  J  df  exp(+i2nfx)  *(v,f)  .  (7b) 

Therefore,  the  region  of  essential  contribution  of  x(v*T)  is 

| v |  <  F,  | t |  <  T,  from  (7b)  and  (7a),  re;  ^ectively. 

The  extents  of  all  four  of  these  two-dimensional  time- 

frequency  representations  are  summarized  in  figure  1.  In  fact, 

2  2 

for  Gaussian  waveform  s(t)  =  a  exp(-ist  /o  ),  the  choices  T  =  4a 
and  F  =  2/(na),  for  example,  give  these  exact  results  in  figure 
1,  at  the  exp(-4)  =  .018  level.  Horizontal  movement  in  this 
figure  is  accomplished  by  means  of  a  Fourier  transform  between 
variables  t  and  \>;  vertical  movement  utilizes  a  Fourier  transform 
relationship  between  x  and  f.  Relations  (6)  and  (7),  along  with 
their  inverse  Fourier  transforms,  constitute  the  totality  of 
these  one-dimensional  transforms. 
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Figure  1.  Extents  of  the  Time-Frequency  Representations 
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GENERALIZED  TIME-FREQUENCY  REPRESENTATIONS 

Since  there  are  four  two-dimensional  domains  of  interest  in 
the  TFRs  depicted  in  figure  1,  it  is  necessary  to  consider  the 
effects  of  weighting  and  smoothing  in  all  of  them. 

TWO-DIMENSIONAL  SMOOTHING  OPERATIONS 

Consider  \>,t  weighting  (or  kernel)  v(v,r)  applied 
multiplicatively  to  CAF  x(v/T>  to  yield  modified  (weighted)  CAF 

X(^>/T  )  =  X(v/T)  v(v,t )  .  (8) 

The  three  equivalent  descriptors  to  weighting  v(v,r),  in  the 
remaining  domains,  are  given  by  Fourier  transform  relations 

V( v , f )  =  J  dr  exp(-i2rtfr)  v(v,r)  ,  (9) 

v(t,t)  =  J  dv  exp(+i2nvt)  v(v,x)  ,  (10) 

V(t,f)  =  J  dr  exp(-i2nfr)  v(t,r)  = 

=  J  dv  exp(+i2nvt)  V(v,f)  = 

=  JJ  dv  dr  exp( +i2nvt-i2nf t )  v(v,x)  .  (11) 


The  last  function,  V(t,f)  in  (11),  will  be  called  the  smoothing 
function,  for  reasons  to  be  seen  below.  The  notational 
convention  adopted  here  is  that  a  Fourier  transform  from  t  to  v 
is  indicated  by  a  tilda,  while  a  Fourier  transform  from  r  to  f  is 
indicated  by  a  capital. 
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GAUSSIAN  EXAMPLE 

Probably  the  simplest  example  of  a  unimodal  two-dimensional 
smoothing  operation  in  ail  four  domains  is  furnished  by  the 
following  Gaussian  example,  where  B  and  D  are  arbitrary: 


v(\>,x)  =  exp(  -Jiv2/B2-nx2/D2 )  ,  (12) 
V( v , f  )  =  D  exp( -nv2/B2-nD2f 2 )  ,  (13) 
v( t,r )  =  B  exp( -nB2t2-nr2/D2 )  ,  (14) 
V(t,f)  =  BD  exp( -nB2t2-nD2f 2 )  .  (15) 


The  effective  areas  of  these  four  two-dimensional  functions,  at 
the  1/e  contour  level  relative  to  each  peak,  are  BD,  B/D,  D/B, 
and  1/ ( BD ) ,  respectively.  It  is  seen  from  (12)  that  B  and  D  are 
the  essential  (positive)  extents  of  weighting  v(\>,t)  in  the  v  and 
t  directions,  respectively.  That  is,  v(B,0)  =  v(0,D)  =  exp(-n) 

=  .043  <<  1  =  v(0,0).  Similarly,  from  (15),  1/B  and  1/D  are  the 
essential  (positive)  extents  of  smoothing  function  V(t,f)  in  the 
t  and  f  dimensions,  respectively.  These  properties  are 
illustrated  in  figure  2,  where  each  contour  depicted  is  at  level 
exp(-Ji)  =  .043,  relative  to  its  peak.  Shortly,  we  will 
generalize  this  smoothing  function  example  to  allow  for  tilts  in 
the  v,r  and  t,f  planes,  thereby  enabling  better  smoothing 
capability  to  be  applied  to  the  WDF ,  without  loss  of  significant 
information . 
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Figure  2.  Two-Dimensional  Smoothing  Functions 
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MODIFIED  TIME-FREQUENCY  REPRESENTATIONS 

The  effects  of  each  of  the  general  smoothing  functions  in 
(8)  -  (11)  on  the  four  two-dimensional  TFRs  (4)  -  (7)  of  the 
previous  section  are  now  investigated;  see  also  [4;  appendix  F], 
The  resultant  generalized  time-frequency  representations  (GTFRs) 
are  indicated  on  the  left-hand  sides  by  bold  type: 

X(»#t )  s  X( »,T)  v(v,t )  ,  (16) 

*(v,f)  s  J  dr  exp( -i2nfr )  x(^/T)  =  *(V/f)  ©  V(v,f)  ,  (17) 

R(t,r)  =  J  dv  exp(+i2nvt)  X(V'T)  =  R ( t , r )  ©  v(t,r)  ,  (18) 

r  tf 

VJ(t,f)  s  J  dr  exp(-i2nfr)  R(t,t)  =  W(t,f)  ©  V(t,f)  .  (19) 

x 

Here,  ©  denotes  convolution  on  x,  with  all  other  variables  held 
fixed;  thus,  for  example,  (17)  is  Jdf'  #(v,f-f')  V(v,f'). 

The  interpretations  of  (16)  -  (19)  are  as  follows:  the  CAF  is 
simply  multiplied  by  weighting  v(v,r);  the  SCF  is  smeared  in 
frequency  f  according  to  V(v,f);  the  TCF  is  smeared  in  time  t 
according  to  v(t,r);  and  the  WDF  is  smeared  in  both  t  and  f 
according  to  smoothing  function  V(t,f).  It  is  this  latter 
two-dimensional  smoothing  (convolution)  operation  in  t,f  space 
that  suppresses  or  eliminates  the  undesired  oscillating 
components  that  are  present  in  the  original  WDF,  at  the  expense, 
of  course,  of  spreading  out  localized  energy  components  of  the 
waveform . 
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The  extents  of  the  GTFRs  are  sketched  in  figure  3;  these 
results  are  based  upon  (16)  -  (19),  in  combination  with  figures  1 
and  2.  Because  x(v/r)  is  the  result  of  multiplication  (16),  its 
extents  in  \>,x  are  the  minima  of  the  two  contributing  functions. 
On  the  other  hand,  the  f  extent  of  *(v,f)  is  increased  by  1/D, 
which  is  the  positive  extent  of  V(v,f)  in  f.  Similarly,  the  t 
extent  of  R(t,x)  is  lengthened  by  1/B,  owing  to  the  smoothing 
action  of  v(t,t).  In  both  of  these  latter  cases,  the  length  of 
the  untransformed  variable  (\>  for  #(v,f)  and  x  for  R(t,x))  is 
unchanged.  Finally,  W(t,f)  is  lengthened  by  1/B  and  1/D  in  the 
t  and  f  dimensions,  respectively,  owing  to  the  double  convolution 
with  smoothing  function  V(t,f). 

Since  the  smoothing  function  V(t,f)  in  figure  2  has 
essentially  reached  zero  by  the  time  |t|  =  1/B  and  |f|  =  1/D,  the 
effective  extents  in  t  and  f  are  approximately  |t|  <  1/(2B)  and 
| £ |  <  1 / ( 2D ) .  That  is,  V(t,f)  is  approximately  1/B  by  1/D  wide 
in  the  t,f  plane,  for  an  effective  area  of  1/(BD);  see  the  line 
under  (15).  If  this  area  1/(BD)  is  .5  or  greater,  then  we  can 
expect  that  smoothed  WDF  W(t,f)  will  be  everywhere  positive 
[4;  ( F-7 )  -  ( F— 19)3. 

On  the  other  hand,  if  effective  area  1/(BD)  is  significantly 
less  than  .5,  then  smoothing  function  V(t,f)  is  rather  impulsive- 
like  and  little  averaging  will  occur  as  a  result  of  double 
convolution  (19).  Thus,  it  appears  that  BD,  at  least  for  the 
simple  Gaussian  example  in  (12)  -  (15)  and  figure  2,  should  be 
chosen  of  the  order  of  3  to  4 .  Then,  the  effective  area  of 
weighting  v(v,x)  in  (12)  and  figure  2  is  BD,  which  is  of  the 


11 


2 


TR  8785 


order  of  3  to  4 .  This  area  is  significantly  smaller  than  the 
effective  extent  of  CAF  \(v,t)  in  figure  1,  which  covers  an  area 
of  the  order  of  FT,  which  is  generally  much  larger  than  1. 
Therefore,  we  can  anticipate  significant  modifications  in  the 
weighted  CAF  x(v'T)/  and,  hence,  in  the  smoothed  WDF  W(t,f),  in 
the  majority  of  the  t,f  plane.  Except  to  say  that  we  expect  that 
B  <  F  and  D  <  T,  there  is  little  quantitative  connection  between 
these  parameters,  in  general. 


TILTED  GAUSSIAN  EXAMPLE 


When  waveform  s(t)  contains  some  linear  frequency  modulation, 
the  simple  Gaussian  smoothing  functions  in  (12)  -  (15)  and  figure 
2  are  inadequate.  The  CAF  and  WDF  of  s(t)  have  contours  in  their 
respective  planes  that  are  similar  to  tilted  ellipses;  see,  for 
example,  [4;  pages  35  -  39].  It  is  then  necessary  to  realize  a 
weighting  function  v(v,r)  and  a  smoothing  function  V(t,f),  which 
also  have  the  capability  of  moving  their  contours  to 
approximately  match  those  of  typical  CAFs  and  WDFs. 

A  very  useful  set  of  smoothing  functions  is  furnished 
by  the  tilted  Gaussian  mountain,  with  B  and  D  arbitrary 
[4;  appendices  F  and  D]: 


2 

2 

1  ■ 

=  exp 

-n 

V 

_2 

+  ^  + 

2r  —  — 

B  D 

IB 

D 

i  . 

V(v,f ) 


D  exp 

-n 

[l-r2)^r-  +  D2 f 2  -  i2r  £  Df 

,  a  H  J 

(21) 
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v(t,r) 


B  exp 

-n 

2  2  ( 
B  t  +  | 

>-r2)?  + 

i2r  Bt  § 

(22) 


V(t,f ) 


ix  fr,2^2  ,  ^2r2  .  ^ 

- =■  B  t  +  D  f  +  2r  Bt  Df 

1-r  ' 


(23) 


For  r  =  0,  these  reduce  to  (12)  -  (15).  Plots  of  weighting 

function  v(v,r)  and  smoothing  function  V(t,f)  are  displayed  in 

figure  4;  the  contours  drawn  are  at  the  exp(-n)  =  .043  level 

relative  to  the  peak  value  of  each  function.  Dimensionless  tilt 

2  h 

parameter  r  is  limited  to  |r|  <  1;  also,  we  define  q  =  (1-r  )  . 

The  smoothing  function  V(t,f)  again  has  essential  extent  2/B 
by  2/D  in  the  t,f  plane;  that  is,  V(t,f)  is  substantially  zero 
for  | t |  >  1/B  or  | f |  >  1/D.  However,  the  effective  area  Atf 
(inside  the  1/e  relative  contour  level)  of  V(t,f)  is  now  q/ ( BD ) , 
which  can  be  considerably  less  than  1/(BD)  for  |r|  near  1,  that 
is,  when  q  <<  1.  Weighting  function  v(v,t)  now  has  essential 
extent  2B/q  by  2D/q  in  the  v,t  plane;  its  effective  area  A^t  is 
BD/q,  which  is  the  reciprocal  of  that  for  smoothing  function 
V(t,f):  A^t  =  1/At^.  Values  of  At^  of  the  order  of  1/3  to  1/4 
are  desired  for  smoothing  purposes;  then,  AVT  ~  3  to  4 . 

Although  effective  area  At^  can  be  considerably  less  than 
1 / ( BD ) ,  the  smearing  caused  by  double  convolution  (19)  still 
leads  to  a  smoothed  WDF  W(t,f)  which  occupies  the  same  region 
indicated  in  figure  3.  The  extents  of  the  four  GTFRs  are  exactly 
the  same  as  figure  3,  except  that  the  limits  on  v  and  r  are  now 


min(F,B/q}  and  min(T,D/q},  respectively;  q 


(1-r2)15 


(24) 
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Figure  4.  Tilted  Smoothing  Functions 
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CHOI -WILLIAMS  KERNEL 

Another  example  of  the  smoothing  operations  in  (8)  -  (11)  is 
furnished  by  [ 5 ] : 

v(v,r)  =  exp{ -v2r 2/cj2  )  ,  a  >  0  ,  (25) 

n^a/ |rj  exp( -n2a2t2/r 2 )  for  t  *  0 
v( t,T )  =  -  ■  ,  (26) 

,  &  ( t )  for  t  =  0, 

n^0/|v|  exp( -n2o2f 2/v2 )  for  v  *  0 
V(v,f)  =  ,  (27) 

,  &  ( f )  for  v  =  0, 

+  00 

V(t,f)  =  2n^a  J  ^  cos(2nvt)  exp( -n2o2f  2/\>2  )  =  (28a) 

0  + 

-f  oo 

=  2n^a  J  ^  cos(2nfr)  exp( -n2o2t2/r2 )  ,  (28b) 

0+ 

provided  that  t  *  0  and  f  *  0.  Integral  (28a)  is  convergent  at 
v  =  0+  only  if  f  *  0  and  is  convergent  at  v  =  +°°  only  if  t  *  0. 
Similarly,  (28b)  converges  at  r  =  0+  only  if  t  *  0  and  converges 
at  t  =  +<°  only  if  f  *  0.  Also,  (28)  yields  V(0,f)  =  00  for  all 
finite  f,  and  V(t,0)  =  ®  for  all  finite  t.  This  smoothing 
function  V(t,f)  in  (28)  has  an  integrable  singularity  all  along 
both  coordinate  axes  since  v(0,0)  =  JJdtdf  V(t,f)  =  1  is  finite. 
Probably  V(t,f)  has  a  logarithmic  singularity  as  tf  +  0.  Letting 
r  =  | t | x  in  (28b),  V(t,f)  is  seen  to  be  a  function  only  of  |tf|. 
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Because  of  these  singularities,  the  actual  numerical 
calculation  of  the  GTFR  W(t,f),  by  means  of  double  convolution 
(19),  appears  very  unattractive;  rather,  the  Fourier  transform  in 
(19)  is  the  recommended  procedure.  The  delta  functions  in  the 
bottom  lines  of  (26)  and  (27)  mean  that 


R( t , 0 )  =  R ( t , 0 )  and  *(0,f)  =  $( 0,f)  .  (29) 

These  results  follow  directly  from  (18)  and  (17),  respectively. 
Therefore,  when  computing  GTFR  R(t,r)  by  means  of  the  Fourier 
transform  in  (18),  the  slice  for  r  =  0  need  not  be  done  at  all, 
but  rather  (29)  should  be  employed.  That  is. 


R( t,t )  = 


r  r  2  2  2 

I  dv  exp(i2nvt)  x(v/T)  exp(-v  r  /o  )  for  r  *  0 


R( t , 0 )  =  |s(t) |  for  r  =  Oj 


.  (30) 


Finally,  GTFR  W(t,f)  is  obtained  by  Fourier  transform  (19). 
Numerous  other  possibilities  for  kernel  v(v,r)  are  listed  in  [1]. 


PRODUCT  KERNELS 

The  weighting  in  (25)  is  an  example  of  a  product  kernel,  that 
is,  the  weighting  takes  the  form 

v(v,r )  =  g ( vr )  ,  0(0)  «  1  .  (30a) 

In  order  that  smoothing  function  V(t,f)  be  real  for  all  t,f,  it 
is  necessary  that  v*(-v,-t)  =  v(v,r)  for  all  v,r,  which  in  turn 
requires  that  g(x)  be  real  for  all  arguments  x.  Now  define 
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G(y)  =  J  dx  exp(-i2rtxy)  g(x)  .  (30b) 

Then  G(-y)  =  G*(y)  for  all  y. 

With  the  help  of  these  functions  and  properties,  we  find  that 
the  rest  of  the  two-dimensional  smoothing  functions  are  given  by 


V(v,f ) 


for  v  *  0 
for  v  =  0, 


v(t,i  ) 


Mt) 


for  r  *  0 
for  r  =  0, 


(30c) 


(  30d ) 


V ( t , f  )  =  2  Re  |  &  exp(i2ntfy)  G (p)  .  (30e) 

0 

This  last  result  shows  that  the  smoothing  function  V(t,f)  for  a 
product  kernel  is  always  a  function  of  the  product  tf,  and  is 
never  a  function  of  t  or  f  separately. 

The  last  integral  on  y  converges  at  y  =  0  if  G(°°)  =  0. 
Alternatively,  it  converges  at  y  =  0  for  G(°°)  *  0  if  tf  *  0.  And 
the  integral  converges  at  y  =  00  if  tf  *  0 . 

On  the  other  hand,  if  tf  =  0,  then  the  last  integral  on  y 
above  is  infinite  if  G { 0 )  *  0;  that  is,  V(t,f)  =  ®  for  tf  =  0, 
which  corresponds  to  both  coordinate  axes  t  =  0  and  f  =  0.  The 
example  in  (25)  is  of  this  nature  and  corresponds  to  the  special 

2  2  L  2  2  2 

case  of  g(x)  =  exp (-xz/a  )  and  G(y)  =  rt^o  exp( -n  y  )  ,  for 
which  G(0)  =  n^a  *  0. 
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DISCRETE-TIME  CONSIDERATIONS 

Up  to  this  point,  it  has  been  assumed  that  s(t)  is  available 
for  continuous  time  t.  Now,  we  address  the  case  where  the  only 
knowledge  of  s(t)  is  through  samples  taken  at  multiples  of  time 
increment  A.  The  proper  treatment  of  these  samples  {s(kA)J,  in 
order  to  obtain  an  unaliased  WDF  W(t,f),  was  determined  in  [2]; 
namely,  it  was  found  necessary  to  take  A  <  1/F,  where  bandwidth  F 
is  specified  in  (3).  Also,  when  an  efficient  FFT  procedure  for 
evaluating  discrete  spectral  values  of  S ( f )  was  employed,  it  was 
found  necessary  to  choose  FFT  size  N  >  2T/A,  where  duration  T  is 
specified  in  (2).  The  following  extension  is  aimed  at  obtaining 
an  unaliased  version  of  smoothed  WDF  W(t,f)  defined  in  (19).  The 
reader  must  be  familiar  with  the  procedures  presented  in  [2]. 


EVALUATION  OF  MODIFIED  CAF  X(v/*) 
As  in  [2;  (69)],  define 


S(f ) 


A  exp(-i2nfAk)  s ( kA )  for  |f|  <  (2A) 
-  k 

0  otherwise 


(31) 


where  the  sum  on  k  is  over  all  nonzero  summand  values.  Then 
since  A  <  1/F,  we  have  S(f)  =  S(f)  for  all  f;  furthermore,  S(f) 
can  be  computed  at  any  f  values  of  interest,  directly  from  the 
available  samples  { s ( kA ) ) .  Therefore,  from  (16),  (7b),  and  (5), 
the  modified  continuous  CAF  is 
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X(v,r)  =  v  ( \> ,  r  )  |  df  exp(i2nfr)  $(v,f)  =  (32a) 

=  v(v,t)  |  df  exp(i2nfr)  S(f+Jjv)  S*(f-*5V)  .  (32b) 


Now,  in  practice,  S(f)  must  be  computed  at  a  discrete  set  of 
points;  in  particular,  when  we  choose  frequency  increment 
Af  =  1/(NA),  where  N  is  arbitrary,  we  obtain 

=  A  Y2  exp( -i2nnk/N)  s(kA)  for  |n|  <  ^  .  (33) 

k  1 

There  is  no  need  to  consider  n  beyond  the  ±N/2  range,  because  the 
argument  f  of  S(f)  then  covers  the  ±1/(2A)  frequency  range,  which 
is  greater  than  the  ±F/2  range  of  S(f)  in  (3).  We  adopt,  as  our 
approximation  to  desired  function  (32),  the  trapezoidal  form 


Xa  < v ' T ) 


»<v,T)  jji  E  exp(i2n^r)  sfgj  +  f)  S*  (jU  -  £ 


for  all  v,r 


(34) 


Now  let  infinite  impulse  train 


Sb(x 


E  -  tti  . 

k 


(35) 


Then,  using  Af  =  1/(NA),  (34)  can  be  expressed  and  developed  as 


Xa<V'T) 


V(V,T  ) 


v( V,T  ) 


V( V,T  ) 


j  df  exp(i2nfi)  $(v,f)  A^  6^  (f)  - 
X( v,r  )  e  X]  6(t  -  jNA) 1  = 

j 

]T  X(v»r  -  jNA )  for  all  v,t  . 

j 


(36) 
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The  sum  on  j  in  (36)  represents  sets  of  aliasing  lobes  spaced  by 
multiples  of  NA  on  the  r  axis.  From  figure  1,  since  the  x  extent 
of  x(v/T)  is  ±T,  the  first  aliasing  lobe  in  (36)  for  j  =  1 
extends  down  to  r  =  NA  -  T.  In  order  that  this  lobe  not  overlap 
the  desired  main  lobe,  j  =0,  we  must  have  T  <  NA  -  T,  or 


N  > 


2 T 
A  7 


(37) 


This  last  constraint  on  Af  is  consistent  with  the  fact  that  the  x 
extents  of  R(t,x)  and  x('J/T)  are  ±T;  see  figure  1  and 
[ 2 ;  page  A-4 ] . 

Equation  (37)  states  that  the  size  of  the  FFT  in  (33)  must  be 
at  least  equal  to  twice  the  number  of  waveform  samples  taken  at 
increment  A  in  duration  T  of  s(t).  When  this  selection  is  made, 
(36)  and  (16)  yield 

Xa(v,r)  =  v(\>,x)  x(^/T)  =  X(^/T)  for  |x|  <  NA/2  ,  all  \>  .  (38) 
That  is,  approximate  GTFR  x=(v/i)/  defined  by  the  sum  in  (34),  is 

a 

equal  to  the  desired  GTFR  within  a  strip  in  the  v,x  plane. 

Now,  in  order  to  convert  (34)  to  a  form  where  we  can  use  the 
spectrum  calculations  (33),  we  limit  v  to  the  values  2n/(NA): 


for  | x |  <  ^|  ,  all  n  . 


(39) 


We  have  dropped  the  subscript  a  on  x(v/T)/  by  virtue  of  (38). 
The  v  increment  in  (39)  is  A^  =  2/(NA),  which  is  less  than  1/T 
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according  to  (37);  this  increment  is  fine  enough  to  track 
variations  of  in  v,  since  the  t  extent  of  the  TFRs  in 

figure  2  is  ±T/2 . 

Finally,  in  order  to  manipulate  (39)  into  an  FFT  form,  we 
restrict  the  r-value  calculations  to  the  set 


2n 
,NA  ' 


mA 


vf^7T,mA)  777-  r  exp( i2njm/N)  S  S*  ^-7777 

vNA  NA  J  '  NAJ  NAJ 


for  | m |  <  ^  ,  all  n 


(40) 


Actually,  since  the  |  \>  |  extent  of  x(v'T)  is  minjF,B/q}  according 
to  (24),  we  only  need  to  consider 


1M 

NA 


<  min{F,B/qJ 


that  is,  | n |  <  ^  min{FA,BA/q}  . 


(41) 


But  since  we  always  have  FA  <  1,  then  |n|  <  N/2  will  always 
suffice.  Thus,  m  and  n  in  (40)  can  be  limited  to  ±N/2.  The  t 
increment  in  (40)  is  A^  =  A  <  1/F,  which  is  consistent  with  the 
fact  that  the  f  extents  of  $(v,f)  and  W(t,f)  are  ±F/2;  see 
figure  1  and  [2;  page  A-4 ] . 

We  have  shown  here  that  if  A  <  1/F  and  N  >  2T/A,  then  an 
unaliased  version  of  GTFR  x(v/r)  i-s  available  and  that  this 
version  can  be  efficiently  computed  by  (40).  These  conditions 
are  the  same  as  those  derived  in  [2;  appendix  D].  The 
multiplication  of  x(v'T)  by  weighting  v(v,t)  in  (16)  or  (32)  to 
obtain  x(v'T)  has  no  effect  on  aliasing  in  the  v,r  plane;  this  is 
an  obvious  result  in  retrospect. 
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EVALUATION  OF  MODIFIED  SCF  *(v,f) 

The  modified  SCF  #(v,f)  is  given  by  (17)  as  the  Fourier 
transform  of  x(v/T)*  Since  x(v'T)  will  only  be  available  at 
increment  A^  =  A,  as  given  by  (40),  we  adopt  as  our  approximation 
the  trapezoidal  form 


*a(v,f)  s  A  XI  exp(-i2nfmA)  ,mb)  = 
m 


=  J  dr  exp(-i2nfr)  X(v/T)  A  &A(t)  = 


f  ,  V 

=  #(v,f)  ©  6 1/A  ( f )  =  XI  *b/f  -  ^  for  all  v,f  . 

m 


(42) 


The  first  aliasing  lobe  for  m  =  1  is  centered  at  f  =  1/A. 

The  f  extent  of  GTFR  #(v,f)  is  If^F  +  1/D),  as  seen  in  figure 
3.  In  order  that  aliasing  be  insignificant  in  (42),  we  must  have 
*iF  +  1/D  <  1/(2A);  that  is,  time  sampling  increment  A  must 
satisfy  the  constraint 


A  < 


1 


(43) 


This  is  tighter  than  the  original  constraint  A  <  1/F,  which  was 
sufficient  for  reconstruction  of  s(t)  and  the  unmodified  TFRs 
such  as  x(v/T)  and  W(t,f).  So  if  we  anticipate  doing  some 
smoothing  of  the  TFRs,  sampling  with  a  time  increment  A 
satisfying  (43)  must  be  undertaken  in  order  to  avoid  aliasing. 

In  this  case,  we  have 

*a(v,f)  =  *(v,f)  for  | f |  <  1/(2A)  ,  all  v  .  (44) 
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As  for  the  actual  evaluation  of  GTFR  *(v,t),  we  use  (42)  and 
(44)  to  get 

=  A  E  exp(  -i2n jm/N)  x(^/mA)  for  |j|  <  f  /  all  v  .  (45) 


Finally,  in  ^rder  to  use  the  available  quantities  in  (40),  we 
restrict  the  calculation  to  the  values 


* 


2n  j 
NA'NA. 


=  A 


E  exp( -i2n jm/N)  X (§§"mA) 
m 

for  | n |  <  |  ,  | j |  <  |  • 


(46) 


This  procedure  in  (46)  yields  unaliased  samples  of  the  GTFR 
*(v,r)  when  (43)  is  satisfied.  It  utilizes  FFT  operations, 
applied  to  the  GTFR  x(v,x),  which  is  available  by  the  FFT 
prescription  in  (40).  The  ranges  of  integers  n  and  j  in  (46)  are 
sufficient  to  cover  the  range  ±1/A  and  ±1/(2A)  in  v  and  f, 
respectively.  But  since  1/A  >  F  +  2/D  by  (43),  the  ranges 
±(F  +  2/D)  and  ±(F/2  +  1/D)  in  v  and  f,  respectively,  are 
adequate  to  fully  cover  the  extent  of  GTFR  *(v,f);  see  figure  3. 
The  increment  A^  =  1/(NA)  in  (46)  is  fine  enough  zo  track  *(v,f) 
in  f,  since  1/(NA)  <  1/(2T)  by  (37),  while  the  r  extent  of  the 
GTFRs  in  figure  3  is  always  less  than  ±T. 
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EVALUATION  OF  MODIFIED  WDF  W(t,f) 

The  modified  WDF  W(t,f)  was  given  by  (19)  as  the  Fourier 
transform  of  R(t,r).  However,  in  analogy  to  the  two  alternatives 
in  (6),  there  is  also  the  form 

W(t,f)  =  J  dv  exp(i2nvt)  *(v,f)  .  (47) 

Since  ♦(  v,  f )  will  only  be  available  at  increment  A^  =  2/'(NA),  as 
given  by  (46),  we  utilize  the  following  trapezoidal  approximation 
to  ( 47 ) : 

Wa(t,£)  =  jjf  E  exp(i2n|2t)  = 

=  [  dv  exp(i2nvt)  *(v,f)  A^  &A  (v)  = 

J  V 

=  W(t,f)  ®  &NA/2(t)  =  E  w(t  -  n^,f)  for  all  t,f  .  (48) 


The  first  aliasing  lobe  for  n  =  1  is  centered  at  t  =  NA/2. 

The  t  extent  of  GTFR  W(t,f)  is  ±(ijT  +  1/B),  as  seen  in  figure 
3.  In  order  that  aliasing  be  insignificant  in  (48),  we  must  have 
‘i'T  +  1/B  <  NA/4;  that  is,  the  FFT  size  N  must  satisfy 


N  > 


2T  _4 
A  BA  * 


(49) 


This  is  more  stringent  than  original  constraint  (37),  which 
sufficed  for  the  unmodified  TFRs.  Again,  an  unaliased  smoothed 
WDF  can  only  be  achieved  if  sampling  increment  A  is  smaller  and 
if  the  FFT  size  N  is  larger,  the  exact  amounts  depending  on  the 
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degree  of  smoothing  desired;  see  figure  2  in  this  regard.  When 
(49)  is  satisfied,  we  have 


f)  = 

=  W(t,f) 

for 

|t|  < 

NA/ 4 

/ 

all  f  . 

nat 

ion 

of  (50), 

(48) 

,  and 

(46) 

now  yields 

w(t 

'NA, 

ii 

M 

exp( 

•  i  2n, 
i2nr7T-t 
NA 

)  •( 

2n 

NA 

-A] 

'naJ 

n 

for 

|t|  < 

NA 

4  ' 

m 

< 

N 

2  • 

(50 


(51) 

Finally,  to  convert  (51)  to  an  FFT,  we  restrict  the  t  values  to 

wf"!A  _i] 

W1  2 ' NA  j 

11 

for  | m |  <  |  ,  | j |  <  |  .  ( 52 ) 


NA 


2  £  exp(i2nnm/N) 


NA '  NA 


Again,  N-point  FFTs  will  realize  the  desired  unaliased  smoothed 
WDF  W(t,f)  at  selected  points  that  cover  its  full  extent,  with 
time  and  frequency  increments  that  track  the  fastest  possible 
variations  of  this  function.  In  particular,  (52)  yields 
A  =  A/2  <  1/ ( 2 F )  and  Af  =  1 / ( NA )  <  1/(2T).  But  since  the 
v  extent  of  *(v,f)  in  figure  3  is  never  larger  chan  ±F,  while  the 
x  extent  of  R(t,r)  is  never  larger  than  ±T,  these  increments  are 
certainly  fine  enough  to  track  the  variations;  also  see 
[2;  appendix  A].  The  ranges  of  integers  m  and  j  in  (52)  cover 
interval  ±NA/4  in  t  and  bandwidth  ±1/(2A)  in  f.  But  since 
”.i/4  >  T/2  +  1/B  and  1/(2A)  >  F/2  +  1/D  according  to  (49)  and 
(43),  respectively,  these  t  and  f  ranges  cover  the  full  extent 
of  smoothed  WDF  W(t,f)  in  figure  3. 


26 


TR  8785 


SUMMARY 

Calculation  of  the  modified  CAF  x(v/T)  can  be  done  without 
changing  the  requirements  on  sampling  increment  A  and  FFT  size  N. 
However,  in  order  to  compute  the  modified  SCF  i(v,f),  the 
sampling  increment  A  must  be  taken  at  a  smaller  value,  in  order 
to  avoid  aliasing  in  frequency  f.  Finally,  in  order  to  compute 
the  modified  WDF  W(t,f),  both  the  sampling  increment  A  must  be 
smaller  and  the  FFT  size  N  must  be  larger,  in  order  to  avoid 
aliasing  in  time  t  and  frequency  f.  If  there  is  interest  in 
calculation  of  the  modified  TCF  R(t,r),  it  can  be  shown  that 
aliasing  will  be  controlled  when  FFT  size  N  satisfies  (49);  the 
constraint  (43)  on  A  need  not  be  met,  insofar  as  calculation  of 
R(t,r)  is  concerned,  although  we  still  need  A  <  1/F.  It  is  only 
when  the  final  transformation  into  the  t,f  plane  is  accomplished 
that  both  constraints  (43)  and  (49)  must  be  met. 

If  integrals  of  products  of  WDFs  or  CAFs  are  of  interest  [6], 
the  rules  on  sampling  rate  and  FFT  size  given  here  should  suffice 
to  get  accurate  numerical  results.  The  aliasing  lobes  have  been 
kept  out  of  the  regions  of  interest,  thereby  minimizing  possible 
interference  effects. 

A  summary  of  the  operations  that  must  be  undertaken  on 
available  time  data  samples  { s ( kA ) }  follows:  compute  the 
spectral  quantities  S  in  (33);  use  these  in  (40)  to  get  samples 
of  the  weighted  CAF  employ  (46)  to  evaluate  the  modified  SCF 

*;  and  use  (52)  to  determine  the  smoothed  WDF  W.  All  of  these 
expressions  use  N-point  FFTs . 
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Since  the  number  of  substantial  samples  of  s(t)  is  T/A 
according  to  (2),  the  FFT  size  N  in  (33)  is  at  least  twice  this 
large;  see  (37)  and  (49).  Thus,  approximately  half  of  the  N 
array  locations  input  to  (33)  will  contain  rather  small 
contributions.  If  s(t)  is  sampled  well  beyond  t  =  ±T/2,  say  for 
|t|  >  T,  these  very  small  values  can  be  "collapsed"  into  the 
available  N  bins  with  no  loss  of  accuracy;  see  [2;  page  5]. 

Candidates  for  weighting  v(v,t)  to  be  used  in  (40)  include 
(12)  or  (20)  or  (25).  The  values  of  parameters  B,  D,  r,  and  a 
will  have  to  be  made  by  inspection  of  CAF  x(v>T)/  which  is  the 
factor  multiplying  v(v,r)  in  (40).  A  check  should  then  be  made 
of  (43)  and  (49)  to  ensure  that  aliasing  is  not  significant. 
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